Consideremos la siguiente integral bidimensional: \[\begin{equation*} I = \int_{1}^{4} \int_{2}^{7} \frac{x^{2}y}{3} \,dx \,dy = \int_{0}^{1} \int_{0}^{1} 5 (5 u_{1} + 2)^{2} (3 u_{2} + 1) \,du_{1} \,du_{2} \end{equation*}\]

A continuación vamos a estimar el valor de \(I\) mediante el método de Montecarlo, tanto de forma directa como usando distintas técnicas de reducción de varianza. Para poder comparar más fácilmente los resultados obtenidos, los recopilaremos en una tabla al final de este documento.

Para cada método estimaremos también su coste en tiempo haciendo uso de las herramientas proporcionadas por el paquete bench (en particular, la función mark analiza el coste en tiempo y en memoria de las expresiones proporcionadas, ejecutando cada una de ellas un cierto número de iteraciones y devolviendo una tabla con distintas medidas, entre ellas la mediana de los tiempos de ejecución de cada iteración). De esta forma, podremos comparar la eficiencia de cada método a la hora de estimar el valor de la integral. Para ello, es fundamental usar la misma unidad de tiempo a la hora de estimar el coste.

unidad_de_tiempo <- "s"

Método directo de Montecarlo

genera_vector_aleatorio <- function() {
  runif(2)
}

g <- function(u) {
  u_1 <- u[[1]]
  u_2 <- u[[2]]
  5 * (5 * u_1 + 2)^2 * (3 * u_2 + 1)
}

n <- 1e4
coste_directo <- bench::mark(
  {
    valores <- replicate(n, {
      u <- genera_vector_aleatorio()
      g(u)
    })
  },
  iterations = 10,
  time_unit = unidad_de_tiempo
)$median

estimacion_directo <- mean(valores)
varianza_directo <- var(valores) / n
eficiencia_directo <- 1 / (varianza_directo * coste_directo)

Muestreo estratificado: asignación proporcional

En primer lugar definimos los estratos y establecemos la forma de generar valores aleatorios dentro de cada estrato. Una manera simple de realizar un muestreo sistemático (es decir, de considerar estratos todos del mismo tamaño) es subdividir el intervalo \((0, 1)\) de forma independiente en cada dimensión.

genera_subintervalos <- function(numero_subintervalos) {
  extremos_derechos <-
    seq_len(numero_subintervalos) / numero_subintervalos
  extremos_izquierdos <-
    c(0, extremos_derechos[-numero_subintervalos])
  data.frame(
    min = extremos_izquierdos,
    max = extremos_derechos
  )
}

# Como en cada dimensión se va a considerar el mismo número de
# subintervalos, basta guardar la información una única vez
numero_subintervalos <- 5
subintervalos <- genera_subintervalos(numero_subintervalos)

# Los estratos vienen dados por todas las combinaciones posibles
estratos <- expand.grid(
  seq_len(numero_subintervalos),
  seq_len(numero_subintervalos)
)
cantidad_estratos <- nrow(estratos)
estratos$probabilidad <- 1 / cantidad_estratos

genera_vector_aleatorio_en_estrato <- function(numero_estrato) {
  numero_estrato_u_1 <- estratos[numero_estrato, 1]
  numero_estrato_u_2 <- estratos[numero_estrato, 2]
  estrato_u_1 <- subintervalos[numero_estrato_u_1, ]
  estrato_u_2 <- subintervalos[numero_estrato_u_2, ]
  c(
    runif(1, min = estrato_u_1$min, max = estrato_u_1$max),
    runif(1, min = estrato_u_2$min, max = estrato_u_2$max)
  )
}

Ahora replicamos el proceso de generar valores en cada estrato, en una cantidad proporcional a su probabilidad, y aplicarles la función g a cada uno de ellos.

n_estratos <- (n * estratos$probabilidad) |> 
  ceiling() |> 
  pmax(2)

coste_estratificado_proporcional <- bench::mark(
  {
    valores <- purrr::map2(
      seq_len(cantidad_estratos),
      n_estratos,
      \(numero_estrato, n_estrato) {
        replicate(n_estrato, {
          u <- genera_vector_aleatorio_en_estrato(numero_estrato)
          g(u)
        })
      }
    )
  },
  iterations = 10,
  time_unit = unidad_de_tiempo
)$median

Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.

estimacion_estratificado_proporcional <- mean(unlist(valores))
varianza_estratificado_proporcional <- valores |> 
  purrr::map_dbl(var) |> 
  weighted.mean(w = estratos$probabilidad) |> 
  magrittr::divide_by(n)
eficiencia_estratificado_proporcional <-
  1 / (varianza_estratificado_proporcional *
    coste_estratificado_proporcional)

Muestreo estratificado: asignación óptima

Consideramos los mismos estratos que antes y, por tanto, la misma forma de generar valores en cada uno de ellos.

Ahora replicamos el proceso de generar en cada estrato una cantidad óptima de valores, determinada mediante un procedimiento en dos etapas, y aplicarles la función g a cada uno de ellos. La estimación del coste en tiempo debe, por tanto, abarcar las dos etapas.

n_tanteo <- 50 * cantidad_estratos

coste_estratificado_optimo <- bench::mark(
  {
    # Estimación de las varianzas de los estratos
    n_estratos <- (n_tanteo * estratos$probabilidad) |> 
      ceiling() |> 
      pmax(2)
    valores <- purrr::map2(
      seq_len(cantidad_estratos),
      n_estratos,
      \(numero_estrato, n_estrato) {
        replicate(n_estrato, {
          u <- genera_vector_aleatorio_en_estrato(numero_estrato)
          g(u)
        })
      }
    )

    # Cantidad óptima de valores en cada estrato
    sigmas <- purrr::map_dbl(valores, sd)
    n_estratos <- 
      (n * estratos$probabilidad * sigmas /
         sum(estratos$probabilidad * sigmas)) |> 
      ceiling() |> 
      pmax(2)

    # Generación de valores en cada estrato
    valores <- purrr::map2(
      seq_len(cantidad_estratos),
      n_estratos,
      \(numero_estrato, n_estrato) {
        replicate(n_estrato, {
          u <- genera_vector_aleatorio_en_estrato(numero_estrato)
          g(u)
        })
      }
    )
  },
  iterations = 10,
  time_unit = unidad_de_tiempo
)$median

Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.

estimacion_estratificado_optimo <- valores |> 
  purrr::map_dbl(mean) |> 
  weighted.mean(w = estratos$probabilidad)
varianza_estratificado_optimo <- valores |> 
  purrr::map_dbl(sd) |> 
  weighted.mean(w = estratos$probabilidad) |> 
  magrittr::raise_to_power(2) |> 
  magrittr::divide_by(n)
eficiencia_estratificado_optimo <-
  1 / (varianza_estratificado_optimo *
    coste_estratificado_optimo)

Muestreo por importancia: primera versión

Consideramos como densidad instrumental el producto de la densidad de una \(\mathrm{Beta}(3, 1)\) para \(u_{1}\) por la densidad de una \(\mathrm{Beta}(2, 1)\) para \(u_{2}\) (tratando de esta forma a las dos coordenadas de manera independiente).

Se puede comprobar gráficamente que la densidad instrumental es «parecida» a \(g(u_{1}, u_{2}) f(u_{1}, u_{2})\).

alfa_instrumental_u_1 <- 3
beta_instrumental_u_1 <- 1
alfa_instrumental_u_2 <- 2
beta_instrumental_u_2 <- 1

densidad_nominal <- function(u) {
  u_1 <- u[[1]]
  u_2 <- u[[2]]
  dunif(u_1) * dunif(u_2)
}

densidad_instrumental <- function(u) {
  u_1 <- u[[1]]
  u_2 <- u[[2]]
  dbeta(u_1,
        shape1 = alfa_instrumental_u_1,
        shape2 = beta_instrumental_u_1) *
    dbeta(u_2,
          shape1 = alfa_instrumental_u_2,
          shape2 = beta_instrumental_u_2)
}

library(plotly)

malla <- data.frame(
  u_1 = seq_len(100) / 100,
  u_2 = seq_len(100) / 100
)

malla |> 
  plot_ly(x = ~u_1, y = ~u_2) |> 
  add_surface(z = outer(
    malla$u_1, malla$u_2,
    Vectorize(function(u_1, u_2) {
      u <- c(u_1, u_2)
      g(u) * densidad_nominal(u)
    })
  )) |> 
  add_surface(z = outer(
    malla$u_1, malla$u_2,
    Vectorize(function(u_1, u_2) {
      u <- c(u_1, u_2)
      densidad_instrumental(u)
    })
  ))

Sin embargo, en la cola de la densidad instrumental la razón de verosimilitud tiende a infinito, sin que esto quede compensado por valores pequeños de \(g\). Esto producirá una varianza infinita al aplicar el método de muestreo por importancia.

malla |> 
plot_ly(x = ~u_1, y = ~u_2) |> 
  add_surface(z = outer(
    malla$u_1, malla$u_2,
    Vectorize(function(u_1, u_2) {
      u <- c(u_1, u_2)
      g(u)
    })
  )) |> 
  add_surface(z = outer(
    malla$u_1, malla$u_2,
    Vectorize(function(u_1, u_2) {
      u <- c(u_1, u_2)
      densidad_nominal(u) / densidad_instrumental(u)
    })
  ))

En primer lugar, establecemos la forma de generar valores aleatorios, que ahora será a partir de la distribución instrumental escogida. Por otra parte, el producto de \(g\) por la razón de verosimilitud (es decir, el cociente entre la función de densidad nominal y la función de densidad instrumental) es el siguiente: \[\begin{align*} g(u_1, u_2) \frac{1}{\frac{u_{1}^{2}}{B(3, 1)} \frac{u_{2}}{B(2, 1)}} &= B(3, 1) B(2, 1) \frac{5 (5 u_{1} + 2)^{2} (3 u_{2} + 1)}{u_{1}^{2} u_{2}} \\ &= \frac{\Gamma(3) \Gamma(1)}{\Gamma(4)} \frac{\Gamma(2) \Gamma(1)}{\Gamma(3)} 5 \Biggl( \frac{5 u_{1} + 2}{u_{1}} \Biggr)^{2} \frac{3 u_{2} + 1}{u_2} \\ &= \frac{5}{6} \Biggl( 5 + \frac{2}{u_{1}} \Biggr)^{2} \Biggl( 3 + \frac{1}{u_2} \Biggr) \end{align*}\]

genera_vector_aleatorio <- function() {
  c(
    rbeta(1,
          shape1 = alfa_instrumental_u_1,
          shape2 = beta_instrumental_u_1),
    rbeta(1,
          shape1 = alfa_instrumental_u_2,
          shape2 = beta_instrumental_u_2)
  )
}

razon_de_verosimilitud <- function(u) {
  densidad_nominal(u) / densidad_instrumental(u)
}

g_por_verosimilitud <- function(u) {
  g(u) * razon_de_verosimilitud(u)
}

# Haciendo uso de la simplificación matemática indicada antes, una
# implementación más eficiente para α_u_1=3, β_u_1=1, α_u_2=2 y β_u_2=1 sería
# g_por_verosimilitud <- function(u) {
#   u_1 <- u[[1]]
#   u_2 <- u[[2]]
#   (5 / 6) * (5 + 2 / u_1)^2 * (3 + 1 / u_2)
# }

Ahora replicamos el proceso de generar valores según la densidad instrumental y aplicarles la función g_por_verosimilitud a cada uno de ellos.

coste_importancia_1 <- bench::mark(
  {
    valores <- replicate(n, {
      u <- genera_vector_aleatorio()
      g_por_verosimilitud(u)
    })
  },
  iterations = 10,
  time_unit = unidad_de_tiempo
)$median

Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.

estimacion_importancia_1 <- mean(valores)
varianza_importancia_1 <- var(valores) / n
eficiencia_importancia_1 <-
  1 / (varianza_importancia_1 * coste_importancia_1)

Muestreo por importancia: segunda versión

Una técnica para conseguir una densidad instrumental con razón de verosimilitud acotada es perturbar ligeramente la densidad nominal. Por ejemplo, puesto que la distribución uniforme se corresponde con una distribución \(\mathrm{Beta}(1, 1)\), podemos considerar como densidad instrumental el producto de dos \(\mathrm{Beta}(1, 1/2)\).

Se puede comprobar gráficamente que la razón de verosimilitud ahora está acotada.

alfa_instrumental_u_1 <- 1
beta_instrumental_u_1 <- 1 / 2
alfa_instrumental_u_2 <- 1
beta_instrumental_u_2 <- 1 / 2

densidad_instrumental <- function(u) {
  u_1 <- u[[1]]
  u_2 <- u[[2]]
  dbeta(u_1,
        shape1 = alfa_instrumental_u_1,
        shape2 = beta_instrumental_u_1) *
    dbeta(u_2,
          shape1 = alfa_instrumental_u_2,
          shape2 = beta_instrumental_u_2)
}

malla |> 
plot_ly(x = ~u_1, y = ~u_2) |> 
  add_surface(z = outer(
    malla$u_1, malla$u_2,
    Vectorize(function(u_1, u_2) {
      u <- c(u_1, u_2)
      densidad_nominal(u) / densidad_instrumental(u)
    })
  ))

El producto de \(g\) por la razón de verosimilitud es ahora el siguiente: \[\begin{align*} g(u_1, u_2) \frac{1}{\frac{(1 - u_{1})^{-1/2}}{B(1, 1/2)} \frac{(1 - u_{2})^{-1/2}}{B(1, 1/2)}} &= B(1, 1/2) B(1, 1/2) \frac{5 (5 u_{1} + 2)^{2} (3 u_{2} + 1)}{(1 - u_{1})^{-1/2} (1 - u_{2})^{-1/2}} \\ &= 20 (5 u_{1} + 2)^{2} (3 u_{2} + 1) \sqrt{(1 - u_{1}) (1 - u_{2})} \end{align*}\]

genera_vector_aleatorio <- function() {
  c(
    rbeta(1,
          shape1 = alfa_instrumental_u_1,
          shape2 = beta_instrumental_u_1),
    rbeta(1,
          shape1 = alfa_instrumental_u_2,
          shape2 = beta_instrumental_u_2)
  )
}

razon_de_verosimilitud <- function(u) {
  densidad_nominal(u) / densidad_instrumental(u)
}

g_por_verosimilitud <- function(u) {
  g(u) * razon_de_verosimilitud(u)
}

# Haciendo uso de la simplificación matemática indicada antes, una
# implementación más eficiente para α_u_1=1, β_u_1=1/2, α_u_2=1 y β_u_2=1/2 sería
# g_por_verosimilitud <- function(u) {
#   u_1 <- u[[1]]
#   u_2 <- u[[2]]
#   20 * (5 * u_1 + 2)^2 * (3 * u_2 + 1) * sqrt((1 - u_1) * (1 - u_2))
# }

Ahora replicamos el proceso de generar valores según la densidad instrumental y aplicarles la función g_por_verosimilitud a cada uno de ellos.

coste_importancia_2 <- bench::mark(
  {
    valores <- replicate(n, {
      u <- genera_vector_aleatorio()
      g_por_verosimilitud(u)
    })
  },
  iterations = 10,
  time_unit = unidad_de_tiempo
)$median

Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.

estimacion_importancia_2 <- mean(valores)
varianza_importancia_2 <- var(valores) / n
eficiencia_importancia_2 <-
  1 / (varianza_importancia_2 * coste_importancia_2)

Muestreo por importancia: tercera versión

Otra técnica, conocida como muestreo por importancia defensivo, consiste en considerar como densidad instrumental una mixtura de la densidad deseada con la densidad nominal. De esta manera tenemos asegurado que la razón de verosimilitud está acotada, ya que \[\begin{equation*} \frac{f_1(x)}{p f_1(x) + (1 - p) f_2(x)} \leq \frac{f_1(x)}{p f_1(x)} \leq \frac{1}{p} \end{equation*}\]

alfa_instrumental_u_1 <- 3
beta_instrumental_u_1 <- 1
alfa_instrumental_u_2 <- 2
beta_instrumental_u_2 <- 1

densidad_instrumental <- function(u) {
  u_1 <- u[[1]]
  u_2 <- u[[2]]
  f_1 <- densidad_nominal(u)
  f_2 <-
    dbeta(u_1,
          shape1 = alfa_instrumental_u_1,
          shape2 = beta_instrumental_u_1) *
    dbeta(u_2,
          shape1 = alfa_instrumental_u_2,
          shape2 = beta_instrumental_u_2)
  0.7 * f_1 + 0.3 * f_2
}

malla |> 
plot_ly(x = ~u_1, y = ~u_2) |> 
  add_surface(z = outer(
    malla$u_1, malla$u_2,
    Vectorize(function(u_1, u_2) {
      u <- c(u_1, u_2)
      densidad_nominal(u) / densidad_instrumental(u)
    })
  ))

En este caso no es posible obtener una expresión simplificada del producto de \(g\) por la razón de verosimilitud.

genera_vector_aleatorio <- function() {
  if (runif(1) < 0.7) {
    runif(2)
  } else {
    c(
      rbeta(1,
            shape1 = alfa_instrumental_u_1,
            shape2 = beta_instrumental_u_1),
      rbeta(1,
            shape1 = alfa_instrumental_u_2,
            shape2 = beta_instrumental_u_2)
    )
  }
}

razon_de_verosimilitud <- function(u) {
  densidad_nominal(u) / densidad_instrumental(u)
}

g_por_verosimilitud <- function(u) {
  g(u) * razon_de_verosimilitud(u)
}

Ahora replicamos el proceso de generar valores según la densidad instrumental y aplicarles la función g_por_verosimilitud a cada uno de ellos.

coste_importancia_3 <- bench::mark(
  {
    valores <- replicate(n, {
      u <- genera_vector_aleatorio()
      g_por_verosimilitud(u)
    })
  },
  iterations = 10,
  time_unit = unidad_de_tiempo
)$median

Finalmente, estimamos el valor de la integral, la varianza de esa estimación y la eficiencia del método.

estimacion_importancia_3 <- mean(valores)
varianza_importancia_3 <- var(valores) / n
eficiencia_importancia_3 <-
  1 / (varianza_importancia_3 * coste_importancia_3)

Tabla comparativa de resultados

knitr::kable(
  data.frame(
    `Método` = c(
      "Directo",
      "Estratificado proporcional",
      "Estratificado óptimo",
      "Importancia: versión 1",
      "Importancia: versión 2",
      "Importancia: versión 3"
    ),
    `Estimación` = c(
      estimacion_directo,
      estimacion_estratificado_proporcional,
      estimacion_estratificado_optimo,
      estimacion_importancia_1,
      estimacion_importancia_2,
      estimacion_importancia_3
    ),
    Varianza = c(
      varianza_directo,
      varianza_estratificado_proporcional,
      varianza_estratificado_optimo,
      varianza_importancia_1,
      varianza_importancia_2,
      varianza_importancia_3
    ),
    Coste = c(
      coste_directo,
      coste_estratificado_proporcional,
      coste_estratificado_optimo,
      coste_importancia_1,
      coste_importancia_2,
      coste_importancia_3
    ),
    Eficiencia = c(
      eficiencia_directo,
      eficiencia_estratificado_proporcional,
      eficiencia_estratificado_optimo,
      eficiencia_importancia_1,
      eficiencia_importancia_2,
      eficiencia_importancia_3
    )
  ),
  digits = 10
)
Método Estimación Varianza Coste Eficiencia
Directo 279.6914 3.9836253 0.03419515 7.3410305
Estratificado proporcional 279.4703 0.1762562 0.88960109 6.3776451
Estratificado óptimo 279.5725 0.1485181 1.03732708 6.4909000
Importancia: versión 1 283.0876 19.0896688 0.15074517 0.3475027
Importancia: versión 2 278.9978 2.8236003 0.17674469 2.0037818
Importancia: versión 3 277.8903 1.0186354 0.20140209 4.8743562